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Abstract 

The detections of atomic hydrogen, heavy atoms and ions surrounding the 
extrasolar giant planet (EGP) HD209458b constrain the composition, tem- 
perature and density profiles in its upper atmosphere. Thus the observations 
provide guidance for models that have so far predicted a range of possible 
conditions. We present the first hydrodynamic escape model for the upper 
atmosphere that includes all of the detected species in order to explain their 
presence at high altitudes, and to further constrain the temperature and 
velocity profiles. This model calculates the stellar heating rates based on 
recent estimates of photoelectron heating efficiencies, and includes the pho- 
tochemistry of heavy atoms and ions in addition to hydrogen and helium. 
The composition at the lower boundary of the escape model is constrained 

* Corresponding author. Fax: +1 (520) 621 4933 
Email address: tommiOlpl . arizona. edu (T. T. Koskinen) 

Preprint submitted to Icarus October 5, 2012 



by a full photochemical model of the lower atmosphere. We confirm that 
molecules dissociate near the 1 //bar level, and find that complex molecular 
chemistry does not need to be included above this level. We also confirm 
that diffusive separation of the detected species does not occur because the 
heavy atoms and ions collide frequently with the rapidly escaping H and 
This means that the abundance of the heavy atoms and ions in the 
thermosphere simply depends on the elemental abundances and ionization 
rates. We show that, as expected, H and O remain mostly neutral up to at 
least 3 whereas both C and Si are mostly ionized at significantly lower 
altitudes. We also explore the temperature and velocity profiles, and find 
that the outflow speed and the temperature gradients depend strongly on 
the assumed heating efficiencies. Our models predict an upper limit of 8,000 
K for the mean (pressure averaged) temperature below 3 Rp, with a typical 
value of 7,000 K based on the average solar XUV flux at 0.047 AU. We use 
these temperature limits and the observations to evaluate the role of stellar 
energy in heating the upper atmosphere. 

Keywords: Extra-solar planets, Aeronomy, Atmospheres, composition, 
Photochemistry 



1 1. Introduction 



The detection of hot atomic hydrogen in the upper atmosphere of HD209458b 



3 (Vidal-Madjar et al. , 2003, 2004) has inspired numerous attempts to model 



4 physical and chemical processes in highly irradiated atmospheres, including 



2 



5 rapid escape as one of the most challenging aspects. Subsequent detection of 



6 heavy atoms and ions ( Vidal-Madjar et al. , 2004 Linsky et al. , 2010) point 



out the need for more complex models that include the chemistry associated 
with these species as well as the collision coupling between them and the ma- 
jor species. Indeed, close-in extrasolar planets offer a natural laboratory to 
constrain the theory of rapid escape, including hydrodynamic escape. This 
is important because aspects of the theory are controversial, and yet rapid 
escape is believed to have played a role in shaping the early evolution of the 



ip,tmospheres in the solar system (e.g., Zahnle and Kasting, 1986 Hunten et 
ahj 1987). Escape may also be a crucial factor in determining atmospheric 



15 conditions and habitability of super-Earths and Earth-like planets around 



16 M dwarfs (e.g.. Tarter et al. , 2007) that may be amenable to spectroscopic 



studies in the near future (e.g., Charbonneau et al. , 2009). 

The basic ideas about the nature of the upper atmospheres around close-in 



19 EGPs were laid out almost as soon as the first planet, 51 Peg b (Mayor et al 



1995 ) , was detected. For instance, Coustenis et al. ( 1998 ) argued that heating 



by the stellar EUV radiation and interaction with the stellar wind leads to 
high temperatures of several thousand Kelvins in the upper atmosphere and 
exosphere of close-in EGPs. They also suggested that the upper atmosphere 
is primarily composed of atoms and ions, and that hydrodynamic escape 
might be able to drag species heavier than H and He into the exosphere. At 



26 the same time, Schneider et al. (1998) argued that material escaping from 



27 the atmospheres of close-in EGPs would form a potentially observable comet- 



3 



28 like tail. When IVidal-Madjar et al.l (120031 [20041) detected the transits of 



HD209458b in the stellar FUV emission lines, they also argued that the planet 
is followed by comet-like tail of escaping hydrogen, and that hydrodynamic 
escape is required to drag oxygen and carbon atoms to the thermosphere. 



The model of Yelle (2004, 2006) was the first attempt to model the aeron- 



omy and escape processes in detail and most of the assumptions in that 
work have been adopted by subsequent investigators. It solved the vertical 
equations of continuity, momentum, and energy for an escaping atmosphere, 
including photochemistry in the ionosphere and transfer of stellar XUV ra- 
diation. Based on a composition of hydrogen and helium, the results demon- 
strated that H2 dissociates in the thermosphere, which at high altitudes is 
dominated by H and H"*". The model also showed that stellar heating leads 
to temperatures of ~ 10,000 K in the upper atmosphere, and predicted an 



energy-limited mass loss rate of 4.7 x 10'' kg s ^ (Yelle, 2006). 



Yelle (2004) argued that conditions beyond ~3 Rp were too complex and 



uncertain to be modeled reliably and therefore chose an upper boundary at 3 
Rp, rather than at infinity, as adopted in early solar wind models. This led to 
a requirement for boundary conditions for the fluid equations at a finite ra- 



46 dius. Yelle (2004) required consistency between fluid and kinetic simulations. 



47 based on the well established fact that kinetic and fluid approaches provide 



consistent results for the escape flux (e.g., Lemaire and Scherer, 1973). This 



led to subsonic velocities of a few km s ^ in his model - although the presence 
of a sonic point at a higher altitude was not ruled out. 
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Many other models for the upper atmospheres of close-in EGPs have 



52 been published (e.g., Lammer et al. , 2003 Lecavelier des Etangs et al. , 2004 



Jaritz et alj [2005| [Tian et al.j [2005| [Erkaev et alj [20071 l^larcia Munozj [20071 



54 l^chneiter et all [20071 Peiiz et al.j [20081 [Holstrom et al.j [20081 |Murray-Clay et 



al. , 2009 Stone and Proga 2009 Guo 2011 Trammell et al. , 2011 ). These in 



elude one-dimensional, two-dimensional, and three-dimensional models that 
make different assumptions regarding heating efficiency, the effect of stellar 
tides, photochemistry, and the escape mechanism. Despite significant dif- 
ferences in the temperature and velocity profiles, almost all of the existing 
models agree that close-in EGPs such as HD209458b are surrounded by an 
extended, hot thermosphere that is undergoing some form of escape. Most of 
the models to date concentrate on the distribution of H and in the upper 



63 atmosphere. Garcia Munoz (2007) developed the only model to address the 



presence of O and C"*" in the thermosphere ( Vidal-Madjar et al. , 2004; Linsky 



et al. , 2010). This model is otherwise similar to Yelle (2004), but it includes 



the photochemistry of heavy ions, atoms, molecules, and molecular ions. It 
also extends to higher altitudes, and includes the effect of substellar tidal 
forces and stellar wind, albeit in an approximate manner. 

Koskinen et al. ( 2007a[ [b) developed a three-dimensional model for the 



thermospheres of EGPs at wide orbits. They pointed out that the atmo- 
spheres of close-in EGPs do not escape hydrodynamically unless they receive 
enough stellar XUV energy to dissociate molecules in the EUV heating layer 
below the exobase. Although their results are limited to the specific case of 
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H2, they can be generalized as follows. The most important molecules H2 
(through the formation of H^), CO, H2O, and CH4 act as strong infrared 
coolants in the thermosphere. High temperatures and rapid escape are only 



possible once these molecules are dissociated. Koskinen et al. (2007b) showed 



that H2 dissociates in the thermosphere of a Jupiter-type planet orbiting a 
Sun-like star within 0.2 AU. Once H2 dissociates, it is reasonable to assume 
that other molecules dissociate too. At this point the pressure scale height 
is enhanced by a factor of ~10 when H becomes the dominant species in the 
thermosphere and temperatures reach 10,000 K. 

It should be noted that a composition of H and H"*" with high temper- 
atures does not guarantee that the atmosphere escapes hydrodynamically. 



85 For instance, Koskinen et al. (2009) showed that hydrodynamic escape is ex- 



tremely unlikely to occur on a planet such as HD17156b because of its high 
mass and eccentric orbit. These types of results have implications on statis- 
tical studies that characterize the escape of planetary atmospheres by relying 



ipn the so-called energy-limited escape (e.g., Watson et al. 1981 Lecavelier 



des Etangs, 2007; Sanz-Forcada et al. , 2010). These studies often include an 



91 efficiency factor in the mass loss rate that is based on the heating efficiency 



of the upper atmosphere (e.g., Lammer et al. , 2009). Unless the atmosphere 
is escaping rapidly, the heating efficiency could be considerably larger than 
the fraction of energy that actually powers escape through adiabatic cool- 
ing. Under diffusion-limited escape or in the Jeans regime the energy-limited 
escape rate is just an upper limit and the true escape rate can be lower. 
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Ideally, the uncertainties in the models can be limited by detailed obser- 
vations of the escaping species. At present, multiple observations are only 
available for HD209458b, and they reveal the presence of H, O, C+, and Si^"*" 



high altitudes in the thermosphere ( Vidal-Madjar et al. 2003 2004; Lin 



sky et al. , 2010). Visible and infrared observations have also revealed the 



ipresence of Na, H2O, CH4, and CO2 in the lower atmosphere (Charbonneau 



et al. , 2002 Knutson et al. 2008 Swain et al. 2009). Taken together, these 



observations are beginning to reveal the composition and thermal structure 
in the atmosphere of HD209458b. The purpose of the current paper is to 
characterize the density profiles of all of the detected species in the ther- 
mosphere, and to explain the presence of the heavy atoms and ions at high 
altitudes in the upper atmosphere. The results can be used to infer some 
basic properties of the atmosphere. 

To this end, we introduce a one-dimensional escape model for the upper 
atmosphere of HD209458b that includes the photochemistry of heavy atoms 
and ions. As pointed out above, previous models agree broadly on the quali- 
tative nature of the thermosphere but the temperature, density, and velocity 



profiles predicted by them differ significantly (see Section 3.1). Some authors 
have argued that the density of H in the thermosphere is not sufficient to 



explain the observed transit depths (see Koskinen et al. , 2010a, for a review 



117 thus lending support to alternative interpretations of the observations such 



118 as the presence of energetic neutral atoms (Holstrom et al. , 2008) or a comet 



119 like tail of hydrogen shaped by radiation pressure (Vidal-Madjar et al. 2003 ). 
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Accurate modeling of the thermosphere is required to enable better judgment 
between different explanations of the observations. 

The differencies between previous models arise from different assumptions 
regarding heating rates and boundary conditions. In addition to modeling 
the density profiles of the detected heavy species, we have improved these 
aspects of the calculations in our work. For instance, the lower boundary 
conditions are constrained by results from a detailed photochemical model 
of the lower atmosphere (Lavvas et al., in preparation) . With regard to the 
upper boundary conditions, we demonstrate that for HD209458b the extrap- 



olated 'outflow' boundary conditions (e.g., Tian et al. , 2005) are consistent 



with recent results from kinetic theory (Volkov et al. , 2011a|[b ) as long as 
the upper boundary is at a sufficiently high altitude - although uncertainties 
regarding the interaction of the atmosphere with the stellar wind may limit 
the validity of both boundary conditions. We highlight the effect of heat- 
ing efficiency and stellar flux on the density and temperature profiles, and 
constrain the likely heating rates by using photoelectron heating efficiencies 



136 based on the results of Cecchi-Pestellini et al. (2009) and our own estimates 



(Section 3.1). As a result we provide a robust qualitative description of the 
density profiles, and constrain the mean temperature and velocity profile 



in the thermosphere. A second paper by Koskinen et al. (2012) (Paper II) 



140 compares our results directly with the observations. 
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141 2. Methods 

142 2.1. Hydrodynamic model 

143 We use a one-dimensional escape model for HD209458b {Rp = 1.32 Rj, 



144 Mp = 0.69 Mj, a = 0.047 AU) that is similar to the models of Yelle (2004) 



145 and Garcia Munoz (2007). Because such models are extensively discussed 



146 in the literature, we include only a brief overview of the model here, with 

147 the emphasis on how it differs from previous work. The model solves the 

148 one-dimensional equations of motion for an escaping atmosphere composed 

149 of several neutral and ionized species: 



dps 1 d 2 
81 pE) 1 9 , 2 „ , 



^Rst 



dp 



n 1 ^ 



dr \ dr J ^ 



(1) 
(2) 

(3) 



150 where ps is the density of species s, f is the vertical velocity, is the diffusive 

151 flux of species s, Rst is the net chemical source term for species s, is a 

152 force term arising from viscous acceleration, E = c^T is the specific internal 

153 energy of the gas, Qr is the specific net radiative heating rate, k is the 

154 coefficient of heat conduction, and $^ is the viscous dissipation functional 

155 (e.g., O'Neill and Chorlton[ 1989). The total density and pressure are given 

156 by p = Yls Ps ^i^d p = UskT, respectively, where electrons contribute to 
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the total pressure. 

We assumed equal temperatures for the neutral species, ions and elec- 
trons, and calculated the electron density at each altitude from the require- 
ment of charge neutrality. The model solves separate continuity equations 
for each species, but treats the atmosphere otherwise as a single fluid. The 
differences in the velocities of the individual species are taken into account by 
including the diffusive flux Fg. We calculated the fluxes by solving simulta- 
neous equations for multiple species based on the diffusion equation given by 



Chapman and Cowling (1970) (equation 18.2,6, p. 344). We also included a 
force term due to the ambipolar electric field given by eE = —{l/ne)dpe/dr, 
where the subscript e refers to electrons, that can be important in highly 
ionized flows. The collision terms account for neutral-neutral, resonant and 
non-resonant ion-neutral, and Coulomb collisions. This method is in princi- 



170 pie similar to those of Yelle (2004) and Garcia Munoz (2007). We verify that 



the single temperature and diffusion approximations are valid for HD209458b 
based on our results in Section 13.2.21 

The model includes heat conduction and terms due to viscosity in both 
the momentum and energy equations. Thus the equations are consistent 
with the level of approximation in the Navier-Stokes (NS) equations. The 
NS equations themselves are a simplification of the 13-moment solution to the 



Boltzmann equation (e.g., Gombosi, 1994) that is valid when the Knudsen 



number Kn = A/L << 1, where A is the mean free path and L is the typ- 
ical length scale for significant changes in density or temperature. Broadly 
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180 speaking, the equations are valid below the exobase, and terms due to heat 

181 conduction and viscosity gain significance as Kn — )• 1. We note that the 

182 exobase on HD209458b is typically located at a very high altitude (see Sec- 

183 tion 



3.1), and viscosity and heat conduction are not particularly important. 



184 We included species such as H, H+, He, He+, C, C+, O, 0+, N, N+, 

185 Si, Si^, Si^^, and electrons in the hydrodynamic model. We also generated 

186 simulations that included Mg, Mg"*", Na, Na"*", K, K"*", S, and S"*", but the 

187 presence of these species did not affect the density profiles of H, O, C+, or Si^"*" 

188 significantly. The model includes photoionization, thermal ionization, and 

189 charge exchange between atoms and ions. The reaction rate coefficients for 

190 these processes are listed in Table [T] Multiply charged ions were included only 

191 if the ionization potential of their parent ion was sufficiently low compared to 

192 the thermal energy and radiation field in the upper atmosphere. We note that 

193 our model also includes impact ionization by thermal electrons. In general, 

194 this can be important for species with low ionization potential such as alkali 



195 metals (e.g., Batygin and Stevenson, 2010), although we find photoionization 



196 to be more significant in the thermosphere (see Section 3.2). 

197 In order to simulate photochemistry in a numerically robust fashion, we 

198 coupled the dynamical model with the ASAD chemistry integrator developed 



199 at the University of Cambridge (Carver et al. , 1997). In most cases we used 

200 the IMPACT integration scheme that is provided by ASAD. We did not 

201 include any molecules in the present simulations, and thus placed the lower 



boundary of the hydrodynamic model at po = 1 /^bar (see Section 2.1.1). 
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Molecular chemistry is not significant in the thermosphere, where our results 



agree qualitatively with Garcia Munoz (2007) despite simpler chemistry (see 



205 Section 3.2). This is an important result because it implies that complex 

206 molecular photochemistry does not need to be included in the models for the 

207 thermosphere. However, the chemistry of molecular ions may be important 

208 on HD209458b below the 0.1 /ibar level and it needs to be studied in greater 

209 detail. 



Table 1: Reaction rate coefficients 



Rate (cm 



Refer 



PI 

P2 

P3 

P4 

P5 

P6 

P7 

Rl 

R2 

R3 

R4 

R5 

R6 

R7 

R8 

R9 

RIO 

Rll 

R12 

R13 

R14 

R15 

R16 

R17 

R18 

R19 

R20 

R21 
R22 
R23 
R24 
R25 
R26 



H 



H 



+ 



He + hi^ -* He + 
O + hi> -» H+ + 
hi/ -> C+ + 
hv ~i N+ + 
hu -> Si+ 4 



C + 
N + 
Si + 
Si+ 
H+ - 
Ho+ 



H + 
4 Ho 



hf 



H + e->H++e + e 



He 

H 4 
H+ 



He+ - 
f He - 



Hc+ 4- 
^ H + 



He 



H 4- He 



4- 



- e 4- e 



O 4- e ^ 0+ - 
C + e-»C++e4-e 
0+ 4- e -> O + h;/ 
C+ + £ ^ C + hu 
C+ + H ^ C 4- H + 
H 
■ He 



C 4- H+ ^ C+ ■ 
C 4- He+ C+ 
0+ 4- H -> O + 
O + H+ -> 0+ 
N 4- e -» N+ 4- e 
N+ 4- e -» N 
Si 4- e ^ Si+ 
Si+ 4- e ^ Si 
Si+ + e -» E. 

Si^+ 4- e ^ Si+ - 
H+ 4- Si ^ H + 
He+ + Si -» He 
C+ + Si ^ C + 
H + Si^+ -» H+ 
H+ 4- Si+ -> H 



H+ 
- H 



f hv 
24- , 



Si+ 
f Si+ 
Si+ 



Si 



2 + 



2.91 X 10" 
1.75 X 10- 



4.0 X 10" 
4.6 X 10- 



■^(300/Te)" 
^(300/Te)° 



I U°-^^ exp(-l7) , U = 2A.6/E^{eV) 



V0.1804-C/ ) 
1.25 X 10- 

1.75 X 10-"(300/T)"-'"exp(-128,000/T) 
3.59 X 10-** ( „ , C/"-^* exp(-C7) , U = 13.6/E^(eV) 
rO-25 ^ 



6.85 X 10 



0.073+C7 , 



I U" 



3.25 X lO-^^CSOO/Te)" '''' 
4.67 X lO-l^CSOO/Te)" "^" 



6.30 X 10-1^(300/T)-1!"3 oxp(-170, OOO/T) 
1.31 X 10-l^(300/T)-''-^12 
2.50 X 10~^^(300/T)-^'^^'' 



J(300/T)- 



' exp(8.6/T) 



7.31 X 10-l''(300/T)-0-23 exp(-226.0/T) 



0.06524-t/ , 

3.46 X io-l^(300/Te)'' '5°'^ 

4.85 X lO-l^CSOO/Te)" "^" 



exp(-!7) , U = 14.5/£;E(eV) 
.2/Be(eV) 



^ 0.6324(7 

1.57 X 'lO-"(300/Te) 
7.41 X 10 

3.30 X 10-" 
2.10 X lO-'' 
2.20 X 10-''(300/T) 



exp{-(7) , !7 = 16.4/Be(eV) 



7.37 X 10' 



-10 



(300/T)- 



0.24 
-0.24 



^Hum mer and Seaton 'l 963| 
JYan et al. 19981 



Verner et al 
Vcrner et al 
Verner et al 
Verner et al 
Verner et al 



1996. 
1996 
1996" 
1996; 
1996' 



^Storey and Hummer 1995 
^Storey and Hummer 1995 



Voronov 


|1997 


Voronov| 


|l997 



Glover and .Jappscn ^2007 
(Voronov 1997 
(Voronov 1997 
IWoodall et al. 20071 
Voodall et al. 200'f) 
Standi et al. 19981 

Standi et al. 1998] 

JGlover and .Jappscn 20071 
^oodall et al. 200'0' 
^^oodall et al. 200f| 
(Voronov 19971 
jAldrova ndTand PcquignQt| |l973[ 



Voronov||l997( 
jAldrova ndi and Pequignotj |T973[ 



19971 



jAldrovandi a nd Pcquiano lj |197 3[ 
JGl ove r and .lappscn 2007[ 
^oodall et al. 2003" 
IWoodall et al. 200?1 

iKingdon and Ferland 1996 I 



210 The upper atmosphere is heated by stellar XUV radiation. We simulated 
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heating and photoionization self-consistently by using the model density pro- 
files and the UV spectrum of the average Sun. The spectrum covers wave- 
lengths between 0.1-3000 A. The XUV spectrum between 0.1-1050 A was 



generated by the SOLAR2000 model (Tobiska et al. , 2000). It includes strong 



215 emission lines separately and weaker lines binned by 50 A. The Lyman a line 



was included with a wavelength spacing of 0.5 A from Lemaire et al. (2005) 



217 and the rest of the spectrum was taken from Woods and Rottman (2002). 



We assumed that most of the Lyman a radiation absorbed by H is reso- 
nantly scattered and does not contribute significantly to the heating of the 
atmosphere. This is because the lifetime of the 2p state of H is only 1.6 ns, 
compared with the typical collision timescale of ~1 s near the temperature 
peak in the thermosphere of HD209458b. 

References for photoabsorption cross sections of the different species are 
included in Table [T] In general, we divided the incident stellar flux by a factor 
of 4 to account for uniform redistribution of energy around the planet. This 
is expected to be approximately valid in the lower thermosphere based on the 



three-dimensional simulations of Koskinen et al. (2010b). In the extended 



upper thermosphere, on the other hand, radiation passes through to the night 
side and leaves only a small region free of direct heating (e.g., see Figure 2 



230 of Koskinen et al. , 2007b). The current model also includes heating due to 



photoabsorption by C, O, N, and metals. This is relatively insignificant - 
although it leads to some additional heating in the lower thermosphere by 
FUV radiation. 
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234 Heating of the thermosphere is mostly driven by photoionization and 

235 the generation of photoelectrons, although direct excitation of atoms and 

236 molecules may also play a role. Photoelectrons excite, ionize, and dissociate 

237 atoms and molecules until they lose enough energy and become thermalized 

238 i.e., share their energy with thermal electrons in Coulomb collisions. Thermal 

239 electrons share their energy with ions and eventually, the neutral atmosphere. 

240 In highly ionized atmospheres such as on HD209458b the photoelectron heat- 



241 ing efficiency can be close to 100 % ( Cecchi-Pestellini et al. , 2009), depending 

242 on the energy of the photoelectrons. We used scaled heating efficiencies that 

243 depend on photoelectron energy to estimate the net heating efficiency in the 



244 atmosphere (Section 3.1 ). 

245 Generally, we refer to two different definitions of heating efficiency in Sec- 



246 tion in order to highlight the effect of heating efficiency on the tempera- 

247 ture and velocity profiles. The net heating efficiency ?7net is defined simply as 

248 the fraction of the absorbed stellar energy that heats the atmosphere. Pho- 

249 toelectron heating efficiency, on the other hand, applies to photoelectrons 

250 with energy Ep = hu — 1^ where Ig is the ionization potential of species 

251 s and hv is the energy of the ionizing photon. The photoelectron heating 

252 efficiency is the fraction of Ep that heats the thermosphere, and it is gen- 

253 erally higher than r^net because it does not account for recombination. The 

254 net heating efficiency is often used to calculate mass loss rates for extrasolar 



255 planets (e.g., Lammer et al. , 2009). Therefore it is important not to confuse 

256 these two definitions of heating efficiency. We included radiative cooling by 
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257 recombination under the assumption that the thermosphere is optically thin 

258 to the emitted photons. This implies that the ionization potential energy Ig 

259 never contributes to heating at any levels. We also considered the influence 

260 of Lyman a cooling by excited H, although this cooling rate is uncertain and 

261 likely to be low for HD209458b. We discuss the effect of different cooling 

262 rates further in Section 13.11 

263 2.1.1. Lower boundary conditions 

264 As stated above, we placed the lower boundary of the hydrodynamic 

265 model at po = 1 yubar and did not include H2 or other molecules in the model. 

266 This decision was motivated by photochemical calculations for HD209458b 

267 (Lavvas et al., in preparation) that we used to constrain the lower bound- 

268 ary condition. The photochemical model was originally developed for the 

but it was recently expanded 

270 to simulate EGP atmospheres. It calculates the chemical composition from 

271 the deep troposphere (1000 bar) up to the thermosphere above the 0.1 nbar 

272 level by solving the coupled continuity equations for all species based on a 

273 database of ~1,500 reaction rate coefficients and 103 photolysis processes. 

274 Forward and reverse rates are included for each reaction with the ratio of 

275 the rate coefficients defined by thermochemical data. Thus, the results are 

276 consistent with thermochemical equilibrium at deep atmospheric levels but 

277 differences develop at higher altitudes due to photolysis, diffusion, and eddy 

278 mixing. At the lower boundary the chemical abundances of the main species 



atmosphere of Titan (Lavvas et al. , 2008a[[b 
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279 (H, C, N, and O) are set to their thermodynamic equihbrium values and, 

280 depending on the vertical temperature profile and their abundances, species 

281 are allowed to condense. 

282 Figure [T] shows the mixing ratios of H2, H, H2O, O, CH4, CO, CO2, and 

283 C as a function of altitude from the photochemical model. In general, the 



284 results are similar to those of Moses et al. (2011). The H2/H transition is 

285 located near 1 //bar. At lower pressures, the mixing ratio of H2 decreases 

286 rapidly with altitude and falls below 0.1 above the 0.1 //bar level. In agree- 



287 ment with Garcia Munoz (2007), the dissociation of H2 is caused by dissocia- 

288 tion of H2O, which leads to the production of OH radicals that attack the H2 

289 molecules. We note that the exact location of the H2/II transition depends 

290 on the temperature profile and, depending on the thermal structure, it could 

291 occur even below the 1 /ibar level. 

292 The major oxygen-bearing molecules, CO and II2O, have roughly equal 

293 abundances from 10 to 10~^ bar. This is in line with thermochemical equi- 



294 i^ibrium at the temperatures and pressures relevant to HD209458b (Lodders 



and Fegley, 2002). H2O and CO are effectively dissociated at pressures lower 



296 than 0.3 and 0.1 //bar, respectively. We note that molecular abundances 

297 could be significant at 0.1-1 /ibar, and technically the results of the hydro- 

298 dynamic calculations are only valid above the 0.1 /ibar level. The presence of 

299 molecules such as II2, H2O, and CO can lead to enhanced UV heating as well 

300 as efficient radiative cooling by Hj]", H2O and CO in the 0.1-1 /ibar region. 

301 The photochemical model also includes the chemistry of silicon. Due 
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Log,Q(Mixing ratio) 

Figure 1: Mixing ratios of key oxygen and carbon-bearing 
species in the dayside atmosphere of HD209458b (Lavvas et al., 
in preparation). 



302 to condensation into forsterite and enstatite (e.g Visscher et al. , 2010), the 

303 abundance of Si in the observable atmospheres of giant planets was thought to 

304 be negligible and thus the photochemistry of silicon in planetary atmospheres 

305 has not been studied before. The photochemical calculations indicate that, 



306 in agreement with thermochemical calculations (Visscher et al. , 2010), SiO 

307 is the dominant silicon-bearing gas. SiO dissociates in the thermosphere at a 

308 similar pressure level as H2O and CO. We note that the detection of Si^"*" in 

309 the thermosphere implies that silicon does not condense in the atmosphere 

310 of HD209458b (Paper II). 

311 In addition to composition, lower boundary conditions are required for 

312 temperature and velocity. We specified Tq and po at the lower boundary, and 
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used them to calculate po from the ideal gas law. The steady state continuity 
equation PqVqt"^ = Fc, where F^. is the flux constant, was used to calculate the 
velocity vq at the lower boundary during each time step. The flux constant is 
solved self-consistently by the model, and it depends largely on the terms in 
the energy equation. In general we assumed that Tq ^ 1,300 K, which is close 
to the effective temperature of the planet. We note that this temperature 
is largely unconstrained. Radiative transfer models for close-in EGPs (e.g.. 



Showman et al. , 2009 , and references therein) do not account for heating by 



321 stellar UV radiation or possible opacity sources created by photochemistry 



322 (e.g., Zahnle et al. , 2009). Therefore these models may not accurately predict 



the temperature at the base of the thermosphere. 

Sing et al. ( 2008a||b ) used observations of the Na D line proflle to constrain 



the temperature proflle in the upper atmosphere. They suggested that Na 
condenses into clouds near the 3 mbar level, and thus predicted a deep min- 
imum in temperature in this region that is required for condensation. The 
detection of Si^"*" implies that condensation of Na in the lower atmosphere is 
unlikely (Paper II), and thus this result is unreliable. Relying on the same 
observations, Vidal-Madjar et al. ( 2011a|[b ) predicted that the temperature 
increases steeply from 1,300 K to 3,500 K near the 1 /ibar level. However, 
their method to retrieve the temperature relies on the density scale height 
of Na to express the optical depth along the line of sight (LOS). This is not 
consistent with the argument that Na is depleted above the 3 mbar level. 
If such a depletion takes place, the density scale height of Na is not the 
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336 same as the scale height of the atmosphere and it cannot be used to retrieve 

337 temperatures. 



2.1.2. Upper boundary conditions 

Previous models of the thermosphere disagree on the details of the density 



and temperature profiles (e.g., Yelle 2004; Tian et al. , 2005 Garcia Munoz 



2007 Murray-Clay et al. , 2009 ). This is partly due to different boundary con- 



342 ditions, although assumptions regarding the heating rates and composition 



are probably more important (see Section 3.1). Unfortunately, the planetary 
wind equations can have an infinite number of both subsonic and supersonic 
solutions. In time-dependent models, the upper boundary conditions in par- 
ticular can determine if the solution is subsonic or supersonic, and they can 



alter the temperature and velocity profiles (e.g., Garcia Munoz, 2007). The 
choice of proper boundary conditions is therefore important. 

Volkov et al7| ( |2011a|ib| ) studied the escape of neutral atmospheres under 



different circumstances by using the kinetic Monte Carlo (DSMC) method. 
Because the fluid equations are a simplification of the kinetic equations at low 
values of Kn, the hydrodynamic model should ideally be consistent with the 
DSMC results both above and below the exobase. Volkov et al. ( 2011a[ [b) 
found that the nature of the solutions depends on the thermal escape pa- 
rameter Xq = GMpm/kToTQ and the Knudsen number Kuq at the lower 
boundary tq of a region where diabatic heating is negligible. They argued 



357 that hydrodynamic escape is possible when Xq < 2-3 (see also Opik, 1963 



19 



Hunten, 1973). When Xq > 3, on the other hand, the sonic point is at such 



a high altitude that the solution is practically subsonic and with Xq> 6 the 
escape rate is similar to the thermal Jeans escape rate. 

The results of the DSMC calculations can be incorporated into hydro- 
dynamic models with appropriate upper boundary conditions. [Volkov et al. 
( 2011a||b ) suggest that the modified Jeans escape rate, which is based on the 
drifting Maxwellian velocity distribution function, is a good approximation 
of the DSMC results in fiuid models, consistent with Yelle| (2004). In order 
to contrast the modified Jeans upper boundary conditions (hereafter, the 
modified Jeans conditions) with other possibilities, we used them and the 
extrapolated upper boundary conditions (hereafter, the 'outfiow' conditions) 



369 adopted by Tian et al. (2005) and Garcia Munoz (2007) in our simulations. 



In general, we placed the upper boundary at 16 Rp. The impact of the 



boundary conditions is discussed in Section 3.1.5 



We formulated the outfiow conditions simply by extrapolating values for 
density, temperature and velocity with a constant slope from below. For the 
modified Jeans conditions, we calculated the effusion velocity Vs at the upper 



jjDoundary separately for each species by using equation (9) from Volkov et 



376 al. (2011b). Using this equation violates the conservation of electric charge 



at the upper boundary because the small mass of the electrons causes their 
velocity to be much larger than the velocity of the protons. In reality charge 
separation is prevented by the generation of an ambipolar electric field that 
ensures that the vertical current is zero at the upper boundary. This elec- 
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381 trie field eauses the ions to escape faster wliile it slows the electrons down. 

382 Effectively this lowers the escape velocity (wesc = ^/iGMjr) of the ions and 

383 increases the escape velocity of the electrons. 

384 In order to incorporate the ambipolar electric field in the modified Jeans 

385 conditions we expressed the Jeans parameters for ions and electrons as: 



GMmi (peQi 



386 where (pe is the ambipolar electric potential, g^^e is the electric charge and 

387 subscripts i and e stand for electrons and ions, respectively. We used these 

388 Jeans parameters to calculate the effusion velocities for the electrons and ions, 

389 and then solved iteratively for 0e by using the condition of zero current i.e., 

390 ne\qe\ve = 'Ylii^iQi'^i- This approach is consistent with kinetic models for the 

391 solar and polar winds (Lemaire and Scherer, 1971a|[b ). Having obtained the 



392 correct effusion velocities for the charged and neutral species, we evaluated 

393 the mass weighted outflow velocity from: 



V = 

P 



394 and used this velocity as an upper boundary condition. The values of tem- 

395 perature and density that are required for this calculation were extrapolated 

396 from below. As the model approaches steady state, the solution approaches 
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397 a value of v at the upper boundary that is consistent with the modified Jeans 

398 velocity. 



399 2.1.3. Numerical methods 

400 We solved the equations of motion on a grid of 400-550 levels with in- 

401 creasing altitude spacing. The radius r„ from the center of the planet at level 



n is thus given by a geometric series (e.g., Garcia Munoz, 2007): 



n-l 



ri 



(7) 



i=l 



where ri = 1.08 Rp, 6zo = 10 km, and fc = 1.014. We solved the equations of 
motions in two parts, separating advection (Eulerian terms) from the other 
(Lagrangian) terms. The Lagrangian solution is performed first, and all vari- 
ables are updated before advection. We used the fiux conservative van Leer 



scheme (e.g., van Leer, 1979) for advection, and the semi-implicit Crank- 



Nicholson scheme (e.g, Jacobson, 1999) to solve for viscosity and conduction 
in the momentum and energy equations, respectively. We employed a time 
step of 1 s in all of our calculations. Despite the sophisticated numerical 
apparatus, the model is still occasionally unstable, particularly in the early 
stages of new simulations. The primary source of the instabilities are pres- 
sure fiuctuations (sound waves) that are not balanced by gravity. We used 



a two-step Shapiro filter (Shapiro, 1970) periodically to remove numerical 



instabilities. We assumed that a steady state has been reached once the fiux 
constant is constant with altitude and the fiux of energy is approximately 



22 



417 conserved. 



418 3. Results 

419 3.1. Temperature and velocity profiles 

420 In this section we constrain the range of mean temperatures and veloci- 

421 ties based on stellar heating in the thermosphere of HD209458b and similar 

422 close-in EGPs. We discuss the general dependency of the temperature and 

423 velocity profiles on the net heating efficiency and stellar flux, and relate this 

424 discussion to new temperature and velocity profiles for HD209458b that are 

425 based on realistic photoelectron heating efficiencies calculated specifically for 

426 close-in EGPs. This discussion is necessary because the temperature and 

427 velocity profiles from previous models of the upper atmosphere differ signif- 

428 icantly, and the differences affect the density profiles of the observed species 



429 (see Section 3.2). As an example, Figure |2] shows the temperature profiles 

430 based on several earlier models. In addition to boundary conditions, the dis- 

431 crepancies evident in this figure arise from different assumptions about the 

432 heating rates. 

433 3.1.1. General dependency 

434 We note that the thermal structure in the upper atmospheres of the giant 

435 planets in the solar system is not very well understood despite modeling and 

436 observations that are far more extensive than those available for extrasolar 



planets (e.g., Yelle and Miller, 2004). It is therefore useful test the reaction 
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2.0x10" 




Figure 2: Some examples of temperature profiles from earlier 
models of the upper atmosphere of HD209458b. The solid line 



is from Figure 1 of Yelle (2004), the dotted line is from the C2 
model in this work (see Section 3.1.2), the dashed line is the 



atomic hydrogen model from Figure 11 of Tian et al. (2005) 



the dashed-dotted hue is the SP model from Figure 3 of [Garcia 
Munoz ( 2007[ ), and the dashed-triple-dotted hue is from Figure 



1 of 



Murray-Clay et al. (2009) 
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438 of the model to different heating rates and profiles. We used our model to 

439 calculate temperature and velocity profiles based on different heating efficien- 

440 cies and stellar fluxes. These profiles are shown in Figure |3| First, we used 



441 the average solar spectrum (Section 2.1 ) and varied the net heating efficiency 

442 r^net from 0.1 to 1. Second, we multiplied the solar spectrum by factors of 

443 2x, lOx, and lOOx, and used T^net = 0.5. The range of enhanced fluxes covers 

444 solar maximum conditions and stars that are more active than the sun. In 

445 each case we assumed that rj^et is constant and does not vary with altitude. 

446 We used planetary parameters of HD209458b and set the upper boundary to 

447 16 Rp with outflow boundary conditions, and the lower boundary to the 1 

448 yubar level with a temperature of 1,300 K. 

449 A net heating efficiency of 50 % is similar to the heating efficiency in 

450 the Jovian thermosphere ( Waite et al. , 1983[ ) , and we may consider this as 

451 a representative case of a typical gas giant (hereafter, the H50 model). The 

452 maximum temperature in the H50 model is 11,500 K and the temperature 

453 peak is located near 1.5 Rp {p = 0.3 nbar). As rj^et varies from 0.1 to 1, the 

454 peak shifts from 1.4 Rp (0.5 nbar) to 1.9 -Rp (0.1 nbar) and the maximum 

455 temperature varies from 10,000 K to 13,200 K. It is interesting to note that 

456 the temperature profile depends strongly on the heating efficiency but the 

457 location of the peak and the maximum temperature depend only weakly on 

458 r^nct- This is because the vertical velocity increases with heating efficiency, 

459 leading to more efficient cooling by faster expansion that controls the peak 

460 temperature while enhanced advection and high altitude heating flatten the 
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Figure 3: Temperature (a) and velocity (b) profiles in the upper 
atmosphere of HD209458b based on different heating efficiencies 
and stellar XUV fluxes. The sohd hues show models based on 
the average solar flux with ?7nct of 0.1, 0.3, 0.5, 0.8, and 1 (from 
bottom to top). The dashed lines show models with ?7net = 0.5 
and stellar flux of 2x, lOx, and lOOx the solar average flux (in 
order of increasing peak temperature). 
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461 temperature gradient above the peak. As a result, the temperature profile is 

462 almost isothermal when ?7net = 1- 

463 It is also interesting that the temperature profiles in the models that are 

464 based on r^net = 0.5 and the solar flux enhanced by factors of 2-100 differ 

465 from models with the average solar flux and a higher heating efficiency. For 

466 instance, one might naively assume that a model with rjnct = 0.5 and 2x the 

467 average solar flux would be similar to a model with the average solar flux 

468 and rjeta = 1- Suprisingly, this is not the case - despite the fact that the 

469 mass loss rates from these models are identical. This is because of the way 

470 radiation penetrates into the atmosphere - doubling the incoming flux is not 

471 the same as doubling the heating rate at each altitude for the same flux. As 

472 the stellar flux increases further, the temperature peak shifts first to higher 

473 altitudes, and then to lower altitudes so that for the lOOx case the peak is 

474 located again near 1.5 Rp. Despite the hugely increased stellar flux, the peak 

475 temperature is only 18,300 K for the lOOx case. This is again because the 

476 enhanced adiabatic and advective cooling driven by faster expansion control 

477 the temperature even in the absence of efficient radiative cooling mechanisms. 



Koskinen et al. (2010a) suggested that the mean temperature of the ther- 



479 mosphere below 3 Rp can be constrained by observations, and used their 



480 ipmpirical model to fit temperatures to the H Lyman a transit data (Vidal- 



Madjar et al. 2003 Ben-Jaffel 2007, 2008 ). A quantity that can be compared 



with their results is the pressure averaged temperature of the hydrodynamic 
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483 model, which is given by: 



/;;r(p)d(inp) 



484 For ?7net ranging from 0.1 to 1, the pressure averaged temperature below 3 Rp 

485 varies from 6,200 K to 7,800 K for the average solar flux. In the H50 model 

486 the pressure averaged temperature is 7,000 K. We note that Tp is a fairly 

487 stable feature of our solutions - in contrast to the details of the temperature 

488 profile and velocities it is relatively insensitive to different assumptions about 

489 the boundary conditions and heating efficiencies. Obviously, Tp depends on 

490 the stellar flux, although it only increases to 9,800 K in the lOOx case. 



491 Koskinen et al. (2010a) inferred a mean temperature of 8,250 K in the 

492 thermosphere of HD209458b with pq = 1 yubar (the M7 model). Taken to- 

493 gether with our results based on solar XUV fluxes, this implies a relatively 

494 high heating efficiency. Alternatively, with T^net = 0.5 it may imply that the 

495 XUV flux of HD209458b is 5-10 times higher than the corresponding solar 

496 flux. This type of an enhancement is not impossible. The activity level of 

497 the star depends on its rotation rate, and the rotation rate of HD209458 may 



498 be twice as fast as the rotation rate of the sun ( Silva- Valio , 2008). However, 

499 the uncertainty of the observed H Lyman a transit depths accommodates 

500 a range of temperatures, and thus we are unable to derive flrm constraints 

501 on the heating rates from the observations. In general, though, the pressure 

502 averaged temperature provides a useful connection between observations and 
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503 temperatures predicted by models that can be exploited to constrain heating 

504 rates. 

505 The effect of changing the heating efficiency on the velocity profile is 
606 quite dramatic. As rj^et ranges from 0.1 to 1 (with the average solar flux), 

507 the velocity at the upper boundary increases from 2.6 km s~^ to 25 km s^^. 

508 However, the velocity does not increase linearly with stellar flux or without 

509 a bound - in the lOOx case the velocity at the upper boundary is only 30 

510 km s^^. An interesting qualitative feature of the solutions is that the sonic 

511 point moves to a lower altitude with increasing heating efficiency or stellar 

512 flux. With r/nct = 0.1 the isothermal sonic point is located above the upper 

513 boundary whereas with ?7net = 1 it is located at 4 Rp. This behavior is re- 

514 lated to the temperature gradient and it is discussed further in Section 3.1.3[ 

515 Basically the sonic point, when it exists, moves further from the planet as 

516 the high altitude heating rate decreases. 

517 It is now clear that assumptions regarding the heating efficiency and ra- 

518 diative transfer have a large impact on the temperature and velocity profiles 

519 and the results from the previous models reflect this fact (see Figure |2|. The 

520 differences between models have implications on the interpretation of obser- 



521 vations. For instance, Vidal-Madjar et al. (2003) and Linsky et al. (2010) 

522 suggested that the UV transit observations probe the velocity structure of 

523 the escaping gas. Obviously, the nature of this velocity structure depends 



524 ipn the properties of the upper atmosphere. On the other hand, |Ben-Jaffel 



and Hosseini (2010) argued that the observations point to a presence of hot 
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526 energetic atoms and ions within the Roche lobe of the planet. We believe 

527 that it is important to properly quantify the role of stellar heating in creat- 

528 ing the hot, escaping material before other options are pursued. This means 

529 that detailed thermal structure calculations that rely on a proper descrip- 

530 tion of photoelectron heating efficiencies are required. Below we discuss a 

531 new approach to modeling the temperature profile in the thermosphere of 

532 HD209458b and its impact on the velocity and density profiles. 

533 3.1.2. Energy balance and temperatures on HD209458b 

534 In the previous section we discussed models where the net heating effi- 

535 ciency ?7net was fixed at a constant value at all altitudes. In this section we 

536 discuss more realistic models of HD209458b that rely on new approximations 

537 of photoelectron heating efficiency and derive an estimate of ?7net based on 

538 these models. Here we also include radiative cooling from recombination and, 



539 in one case, H Lyman a emissions by excited H (Murray-Clay et al. , 2009). 

540 Our aim was to calculate the most likely range of temperatures in the ther- 

541 mosphere of IID209458b based on average solar fiuxes. Figure |4] shows the 

542 temperature and velocity profiles at 1-5 Rp based on different approxima- 

543 tions (see Table [2] for the input parameters). Model CI assumes a constant 

544 photoelectron heating efficiency of 93 % at all altitudes and photoelectron 

545 energies. This heating efficiency is appropriate for photoelectrons created by 



546 ^0 eV photons at an electron mixing ratio of a;e = 0.1 ( Cecchi-Pestellini et 



547 al.l 2009). Model C2 is otherwise similar to CI but the heating efficiency 
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548 varies with photoelectron energy and altitude (see below). Models C3 and 

549 C4 are also based on CI, but C3 includes the substellar tidal forces in the 



550 equations of motion (e.g., Garcia Munoz, 2007) and C4 includes Lyman a 



551 cooling. All of these models are based on the outflow boundary conditions 

552 for temperature, velocity, and density. 

Table 2: Model input parameters and results 



ModeP 




'/net 


M (10^ kg s-i) 


T, (K)« 


CI 


16 E 


0.56 C 


5.6 


7250 


C2 


16 E 


0.44 V 


4.0 


7200 


C3 


16 E,T 


0.57 C 


6.4 


6450 


C4 


16 E 


0.46 C 


4.5 


7110 


C5 


36 J 


0.56 C 


5.6 


7290 


C6 


16 J 


0.45 V 


3.9 


7310 


s assume po = 10 ^ 


bar and To 


= 1,300 K. 





- Outflow conditions, J - Modified Jeans conditions, T - Substellar tide. 



'^Net heating efficiency (see Section 2.1 ) i.e., the ratio of the net heating flux at ah 
wavelengths to the unattenuated stellar flux (0.45 W m~'^) at wavelengths shorter 
than 912 A. 

'^C - Constant photoelectron heating efficiency, V - Varying photoelectron heating 
efficiency (see text). 

^Pressure averaged temperature below 3 Rp. 



553 Cecchi-Pestellini et al. (2009) also estimated the heating efficiencies for 

554 photoelectrons released by photons of energy > 50 eV at different values 

555 of the electron mixing ratio Xg. We used their heating efficiencies for Xe = 0.1 

556 in the C2 model. They parameterized their results in terms of the vertical 

557 column density A^h of H. We fitted the heating efficiency as a function of Ah 

558 for 50 eV photons with a regular transmission function, and modified this 



31 



559 function accordingly for different cutoff altitudes and heating efficiencies of 



560 photons with different energies [see Figures 3 and 4 of Cecchi-Pestellini et al. 



561 (2009)]. We note that Xg ~ 0.1 near the temperature peak of our models and 

562 thus the results are appropriate for our purposes. However, they are only 

563 applicable to photons with E > 50 eV. We used simple scaling to estimate 

564 the heating efficiencies for low energy photons with E < 50 eV. 

565 As A'h increases, the heating efficiency for 50 eV photons saturates at 93 

566 %. We assumed that the saturation heating efficiency for low energy photons 

567 is also 93 %. In reality, this heating efficiency may be closer to 100 % but 

568 the difference is small. In order to estimate the altitude dependence of the 

569 heating efficiency, we note that the rate of energy deposition by Coulomb 

570 collisions between photoelectrons of energy Ep and thermal electrons with a 

571 temperature T can be estimated from: 



dFE 

dr 



L{Ep,e)%ene [eV cm^^^ s"^] (9) 



572 where Fe is the flux of energy, is the photoelectron flux (cm ^ s Ue 

573 is the density of thermal electrons (cm~'^) and 



, 3.37x10-12 / E„-Ee V-^^ 21 



with Ee = 8.618 x 10 ^Te is the stopping power (Swartz et al. , 1971 ^ 
^Due to a historical precedent, the units here are in cgs. 
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Figure 4: Temperature (a) and velocity (b) profiles in the upper 
atmosphere of HD209458b based on different models (see Table[2] 
for the input parameters). 
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675 Assuming that all of the energy is deposited by electrons that are thermalized 

576 within a path element dr, we can estimate the e-folding length scale for 

577 thermalization of photoelectrons with different energies as follows: 



578 We calculated Ape for different photoelectrons based on the CI model, 

579 and compared the result with the vertical length scale H of the atmosphere. 

580 The latter is either the scale height or Rp, depending on which is shorter. 

581 When Ape/H > 0.005-0.01 we assumed that the heating efficiency decreases 

582 with altitude according to the transmission function for 50 eV photons. The 

583 limiting value was chosen to obtain a rough agreement with the results of 



584 Cecchi-Pestellini et al. (2009) for 50 eV photons, and it implies that the 

585 heating efficiency approaches zero when Ap^,/ H > 0.1. We parameterized the 

586 result in terms of the column density of H based on the density profiles of the 



587 



CI model, and connected our results for low energy photoelectrons smoothly 



with the results of Cecchi-Pestellini et al. (2009) for photons with E > 50 eV. 



589 We then generated the C2 model from the CI model with the new heating 

590 efficiencies. Figure |5] shows the resulting heating efficiencies for 20, 30, 48, 

591 and 100 eV photons. 

592 Figure [6] shows the volume heating rate due to EUV photons of different 

593 energies as a function of pressure based on the C2 model. The maximum 

594 temperature of 12,000 K is reached near 1.5 Rp {p = 0.6 nbar). This re- 
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Figure 5: Heating efficiencies for pfiotons of different energies 
(see text). 

gion is heated mainly by EUV photons with wavelengths between 200 and 
900 A {E = 14-62 eV). The saturation heating efficiency of 93 % for these 
photons is higher than the corresponding heating efficiency in the Jovian 



thermosphere (Waite et al. , 1983). This is because of strong ionization that 



leads to frequent Coulomb collisions between photoelectrons and thermal 
electrons. Radiation with wavelengths shorter than 300 A (E > AO eV) or 
longer than 912 A (13.6 eV) penetrates past the temperature peak to the 
lower atmosphere. The heating efficiency for photons with E > 25 eV ap- 
proaches zero at high altitudes where heating is mostly due to low energy 
EUV photons. The net heating efficiency for the C2 model is ?7net = 0.44 
(Table [2]), which is close to the H50 model. The location of the peak and 
maximum temperature in the C2 model agree with the H50 model, but the 
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Figure 6: Volume heating rate as a function of pressure in the C2 
model due to the absorption of stellar XUV radiation between 
1 and 1000 A in 200 A bins. 

607 temperature at higher altitudes in the C2 model decreases much more rapidly 

608 with altitude than in the H50 model. 

609 Figure [7] shows the terms in the energy equation based on the C2 model. 

610 In hue with previous studies, stellar heating is mainly balanced by adiabatic 

611 cooling. Advection cools the atmosphere at low altitudes below the temper- 

612 ature peak, whereas at higher altitudes it acts as a heating mechanism. In 

613 fact, above 2 Rp the adiabatic cooling rate is higher than the stellar heating 

614 rate because thermal energy is transported to high altitudes by advection 

615 from the temperature peak. The radiative cooling term that is centered near 

616 1.3 Rp arises from recombination following thermal ionization. Recombina- 

617 tion following photoionization is included implicitly in the model and the rate 
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Figure 7: Volume heating rate based on the C2 model (absolute 
values). Advection acts as a cooling mechanism below 1.5 Rp 
and a heating mechanism above this level. 

618 is not included in the output. Conduction is not significant at any altitude 

619 in the model. We note that the rates displayed in Figure [7] balance to high 

620 accuracy, thus implying that the simulation has reached steady state. 

621 The differences between the CI and C2 models are subtle. The peak tem- 

622 peratures are similar, and the temperature profiles generally coincide below 

623 3 Rp. Above this radius the temperature in the C2 model decreases more 

624 rapidly with altitude than in the CI model and subsequently the sonic point 

625 moves to higher altitudes above the model domain. The results indicate that 

626 the assumption of a constant photoelectron heating efficiency is appropriate 

627 below 3 Rp whereas at higher altitudes it changes the nature of the solu- 

628 tion. This should not be confused with the assumption of a constant rj^et, 
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which leads to a different temperature profile when compared with either CI 
and C2 (see Figure |3]). In general, the maximum and mean temperatures 
in models C1-C4 are relatively similar. Thus we conclude that the mean 
temperature in the thermosphere of HD209458b is approximately 7,000 K 
and the maximum temperature is 10,000-12,000 K. 

The substellar tide is included in the C3 model. We included it mainly 



ifo compare our results with previous models (Garcia Munoz, 2007; Penz et 



636 al.l 2008 Murray-Clay et al. 2009). The substellar tide is not a particu 



larly good representation of the stellar tide in a globally averaged sense. In 
reality, including tides in the models is much more complicated than sim- 



639 ply considering the substellar tide (e.g., Trammell et al. , 2011). Compared 



to the CI model, the maximum temperature in the C3 model is cooler by 
~1,000 K and at high altitudes the C3 model is cooler by 1,000-2,000 K. 
The velocity is faster and hence adiabatic cooling is also more efficient. The 



643 substellar tide drives supersonic escape (see also, Penz et al. , 2008 ) and the 



644 sonic point in the C3 model is at a much lower altitude than in the CI model 



645 (see Section 3.1.3). However, it is not clear how the sonic point behaves as a 



function of latitude and longitude. Given that the tide is also likely to induce 
horizontal flows, it cannot be included accurately in ID models. 



Murray-Clay et al. (2009) argued that radiative cooling due to the emis- 



sion of Lyman a photons by excited H is important on close-in EGPs. The 
photons are emitted when the 2p level of H, which is populated by collisions 
with thermal electrons and other species, decays radiatively. We included 
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652 this cooling mechanism in the C4 model by using the rate coefficient from 



Glover and Jappsen (2007) that includes a temperature-dependent correction 



654 to the rate coefficient given by Black (1981). We also included an additional 



655 correction factor of 0.1 based on detailed level population and radiative trans- 



656 fer calculations by Menager et al. (2011). The effect of Lyman a cooling is 



657 largest near the temperature peak where the C4 model is 1,500 K cooler than 

658 the CI model, but generally the difference is not large. We note that the 

659 H Lyman a cooling rate here cannot be generalized as such to other EGPs 

660 because the level populations and opacities depend on the thermal structure 

661 and composition of the atmosphere. 



662 3.1.3. Critical points 

663 As we have pointed out, the location of the sonic point depends on the 

664 energy equation through the temperature proffie. Here we show that the use 

665 of the isothermal approximation in estimating the location of the sonic point 

666 can lead to significant errors unless the atmosphere really is isothermal. The 

667 inviscid continuity and momentum equations can be combined to give an 



668 expression for the critical point of a steady-state solution (Parker, 1965): 



1 dc2 2c2 



(11) 



669 where ^ = r /tq, c = ^JkT /m is the isothermal speed of sound, W = GMp/r^, 

670 and m is the mean atomic weight. It is often assumed that the vertical 

671 velocity at the critical point is given by v"^ = c^(^c) so that the critical point 
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672 coincides with the isothermal sonic point (Parker, 1958). However, Parker 



673 (1965) suggested that subsonic solutions are also possible if the density at 

674 the base of the flow exceeds a critical value determined from the energy 

675 equation. In fact, he argued that conduction at the base of the corona may 

676 not be sufficient to drive a supersonic solar wind. This led him to suggest 

677 that supersonic expansion is possible only if there is significant heating of 



678 



the corona over large distances above the base. 



679 For an isothermal atmosphere with a temperature Tq, equation (11) re 



680 duces to the famous result for the altitude of the sonic point (Parker, 1958): 



2 



W 

682 where W'^/c^ is the thermal escape parameter Xq at the lower boundary 

683 Tq. The isothermal sonic point in the CI model is located at 7.2 Rp where 

684 c(^c) = 7.2 km s~^. The volume averaged temperature of the CI model below 

685 this point is approximately 7,100 K. Assuming that ro = Rp, Tq = 7,100 



686 K, and m = nin, Xq ~ 16 and equation (12) yields C,c ~ 8. In this case 

687 the analytic result agrees fairly well with the hydrodynamic model if one 

688 accounts for partial ionization of the atmosphere by assuming that the mean 

689 atomic weighlj^is m = 0.9 rriH- 

690 On the other hand, the isothermal sonic point in the C2 model is at 15.4 



^The mean atomic weight can be less than 1 because electrons contribute to the number 
density but not significantly to the mass density. 
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691 Rp where c(^c) =4.1 km s^^. This is because the temperature gradient 

692 of the model is steeper than the corresponding gradient in the CI modeL 

693 The volume averaged mean temperature below 15 Rp in the C2 model is 



694 3,900 K. With this temperature and m = rrin, equation (12) predicts a 

695 sonic point at 14.6 Rp. However, at 15 Rp the atmosphere is mostly ionized 



696 and m = 0.6 mn- With this value, the sonic point from equation (12) 

697 would be located at 8.8 -Rp. These examples show that there are significant 



698 caveats to using equation (12) to estimate the altitude of the sonic point on 

699 close-in EGPs without accurate knowledge of the temperature and density 

700 profiles. A variety of outcomes are possible and it is difficult to develop 

701 a consistent criteria for choosing values of T and m that would produce 

702 satisfactory results. 

703 Another problem is that the atmosphere is not isothermal. In fact, the 

704 temperature gradient above the heating peak in models C1-C4 (Table [2]) is 

705 relatively steep, and in some cases it approaches the static adiabatic gradient 



706 (T oc r~^) as defined by Chamberlain ( 1961 ). Assuming that the temperature 

707 profile can be fitted with = Cq/^^ above the heating peak, the estimated 

708 values of (3 for the CI and C2 models are 0.7 and 0.9, respectively. We note 

709 that the velocity in the CI model exceeds the isentropic speed of sound (c^ = 



710 ^J'^kT /m where 7 = 5/3) at 9.8 Rp where = 8.7 km s^^. This altitude 

711 is significantly higher than the altitude of the isothermal sonic point. The 

712 velocity in the C2 model does not exceed the isentropic speed of sound below 

713 the upper boundary of 16 Rp. Thus the temperature profile has a significant 
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714 impact on the nature of the solution and the escape mechanism. This means 

715 that estimating the altitude of the sonic point without observations and 

716 detailed models for guidance is almost certain to produce misleading results. 

717 Past models for the upper atmosphere of HD209458b have predicted a 



718 variety of altitudes for the sonic point. On the other hand, Yelle (2004) 



719 pointed out that stellar heating in the thermosphere is mostly balanced by 



720 adiabatic cooling and our calculations confirm this. Parker (1965) argued 

721 that the critical point stretches to infinity when /3 — )■ 1 i.e., as the temperature 

722 gradient is close to adiabatic. Based on this, we should perhaps expect that 

723 the sonic point on close-in EGPs is located at a fairly high altitude. This 

724 is confirmed by our hydrodynamic simulations. In all of our models except 

725 for one, the sonic point is located significantly above 5 Rp. The exception 

726 is the C3 model, which includes the substellar tide. The isentropic sonic 

727 point in this model is located at 3.9 Rp where = 8.2 km s~^. This is 

728 because the substellar tide leads to a lower effective value of the potential 

729 W . However, the tidal potential depends on latitude and longitude, and the 

730 substellar results cannot be generalized globally. 

731 3.1.4- Mass loss rates 

732 Here we evaluate the mass loss rates based on our models. We define the 

733 mass loss rate simply as: 

M = ATrr'^pv. (13) 
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734 We note that the solar spectrum that we used in this study contains the total 

735 flux of 4 X 10^^ W at wavelengths shorter then 912 A (the ionization 

736 limit of H) when normalized to 1 AU. This value is close to the average solar 



737 fiux of 3.9 X 10^ W m" at the same wavelengths (e.g., Ribas et al. , 2005). 

738 In order to simulate a global average, we divided the fiux by a factor of 4 

739 in the model. This means that the incident flux on HD209458b at 0.047 AU 

740 with wavelengths shorter than 912 A in our model is 0.45 W m~^. The net 

741 heating efficiencies given in Table [2] are based on this value. 

742 Considering first the models with constant r/net ranging from 0.1 to 1 (see 



743 Section 3.1.1), the mass loss rate varies almost linearly with r/not from 10^ 

744 kg s~^ and 10® kg while the pressure averaged temperature below 3 Rp 

745 changes only by 1,500 K. This is because in a hydrodynamic model such 

746 as ours the net energy has nowhere else to go but adiabatic expansion and 

747 cooling, and thus escape is energy-limited. The bulk of the energy is absorbed 

748 below 3 Rp, and the mass loss rate is largely set by radiative transfer in this 

749 region. The mass loss rate for HD209458b predicted by the C2 model is 

750 4.1 X 10^ kg s^^ (r^net = 0.44). The C3 model has the highest mass loss rate, 

751 although this rate is only higher by a factor of 1.13 than the mass loss rate 

752 in the CI model. Thus we predict a mass loss rate of 4-6 x 10"^ kg s~^ from 

753 HD209458b based on the average solar flux at 0.047 AU. 



Garcia Munoz (2007) demonstrated that the mass loss rate is insensitive 



755 to the upper boundary conditions even when they have a large impact on 

756 the temperature and velocity profiles, particularly at high altitudes. Indeed, 
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complex hydrodynamic models are not required to calculate mass loss rates 
under energy- limited escape as long as reasonable estimates of the net heating 
efficiency are available. It is also important to note that the current estimates 



of mass loss rates based on the observations (e.g., Vidal-Madjar et al. , 2003) 



are not direct measurements. Instead, they are all based on different models. 
However, models that predict the same mass loss rate can predict different 
transit depths and models predicting different mass loss rates can match the 
observations equally well. Thus the models should not be judged on how well 
they agree with published mass loss rates but rather on how well they agree 
with the observed density profiles or transit depths. Hydrodynamic models 
with realistic heating rates are required to match the observations, and the 
mass loss rate then follows. 

The globally averaged mass loss rate of about 4-6 x 10^ kg s~^ from 



HD209458b agrees well with similar estimates calculated by Yelle (2004 



2006) and Garcia Munoz (2007), respectively, but it is significantly larger 



than the value calculated by Murray-Clay et al. (2009). These authors re 



port a mass loss rate of 3.3 x 10'^ kg s~^ based on the substellar atmosphere. 
When multiplied by 1/4 this corresponds to a global average rate of about 
8.3 X 10^ kg s~^. However, the substellar mass loss rate is also enhanced by 
tides, so a comparable global average taking this into account would be even 
less than 8.3 x 10^ kg s~^, which is already roughly a factor of 6 smaller than 
our calculations. 



779 The Murray-Clay et al. (2009) models differ in many respects from the 
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models described here including the treatment of boundary conditions and 
radiative cooling, the numerical approach, and the adoption of a gray ap- 
proximation for stellar energy deposition. In order to explore the reason for 
the disagreement in escape rates, we have modified our model to implement 



the gray assumption by using the approach described in |Murray-Clay et aL 
(2009) (see Section 3.2). Specifically, we adopted an incident fiux[^ of 0.45 
W m~^ and a mean photon energy of 20 eV. The mass loss rate based on the 
substellar atmosphere for this model is 2.8 x 10^ kg s^^, in good agreement 



788 with the Murray-Clay et al. (2009) results. Thus, the difference between the 



Murray-Clay et al. (2009) models and the others is due to the gray assump- 



tion, and the fact that they estimated the incident flux on HD209458b based 
on the solar flux integrated between 13 eV and 40 eV. This energy range 
contains only about 25 % of the total solar flux at energies higher than 13.6 
eV. 



794 Although not discussed by Murray-Clay et al. (2009), the restricted en- 



ergy range is likely an attempt to account for the fact that the absorption 
cross section decreases with energy implying that photons of sufficiently high 
energy will be absorbed too deep in the atmosphere to affect escape or the 
thermal structure, or composition of the thermosphere. Whether this is true, 
however, depends on the composition and temperature of the atmosphere. 
The gray assumption also fails to include the fact that the net heating effi- 



■^By chance the incident flux is equal to the mean solar flux divided by 4 that we used 
as a 'globally averaged' value. Here, however, it is taken to be the substellar value. 
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801 ciency increases with higher photon energy. These difficulties highhght the 

802 basic problem with a gray model, that the results can only be accepted with 

803 confidence if verified by a more sophisticated calculation or direct observa- 

804 tions. 

805 3.1.5. Constraints from kinetic theory 

806 Hydrodynamic models should be consistent with kinetic theory of rar- 

807 efied media even if the modeled region is below the exobase. This is because 

808 the atmosphere is escaping to space, and the density decreases with altitude, 

809 falling below the fiuid regime at some altitude above the exobase. There- 

810 fore the conditions in the exosphere affect the flow solutions even below the 

811 exobase. Inappropriate use of the hydrodynamic equations can lead one to 

812 overestimate the flow velocity and mass loss rate, and these errors also affect 

813 the temperature and density proflles. Thus it is important to demonstrate 

814 that the hydrodynamic solutions agree with constraints from kinetic theory 

815 (e.g., [Volkov et al.[ |2011a||b| . 



816 As an example, we calculated Kn^ and Xq (see Section 2.1.2 ) based on the 

817 CI and C2 models. The Knudsen number Kn depends on the mean collision 

818 frequency, and it is much smaller than unity at all altitudes below 16 Rp. 



819 -Thus the exobase is located above the model domain (see also Murray-Clay 



et al. , 2009). Calculating values for Xq is complicated by the broad stellar 



821 heating proflle. We consider the region where stellar heating is negligible to 
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822 be where the flux of energy 



E^ = F^\ c,T + - j - ^r'^ (14) 



823 is approximately constant. This criteria is consistent with the equations of 

824 motion, and it means that tq that should be used to calculate Xq is above 

825 the upper boundary of our model because signiflcant stellar heating persists 

826 at all altitudes. Thus we evaluated values of X near the upper boundary for 

827 guidance. We also calculated the values with both the mass of the proton 

828 (Xh) and the mean atomic weight {X). 

829 In the CI model, X^ decreases with altitude, and above 11.4 Rp it has 

830 values of less than 3. The mean atomic weight near the upper boundary is 

831 ~0.6 amu, and thus the general value of X < 2 above 11.1 Rp. The sonic 

832 point in the CI model is below 11 Rp, and it is in a region where stellar 

833 heating is signiflcant. In the C2 model, both X and Xh are greater than 3 

834 at all altitudes below 16 Rp. In fact, X increases with altitude above 9 Rp 

835 because the temperature gradient parameter exceeds unity. Thus the values 

836 X in the CI and C2 models are consistent with the difference in altitude 



837 between the sonic points in these models (see Section 3.1.3). Indeed, our 

838 results show, in line with Parker's original ideas about the solar wind, that 

839 supersonic escape is possible if there is signiflcant heating of the atmosphere 

840 over large distances above the temperature peak. Such heating flattens the 

841 temperature gradient and brings the sonic point closer to the planet. 
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842 We note that there are some caveats to applying the simple criteria based 

843 on Kuq and Xq to close-in extrasolar planets. The upper atmospheres of 



844 ifhese planets are strongly ionized, and the DSMC simulations of Volkov et 



845 |aL] ( 2011a||b ) apply only to neutral atmospheres. Partly due to ionization, 

846 the collision frequencies in the thermospheres of close-in planets are also 

847 high. Further, the atmospheres are affected by a broad stellar heating pro- 

848 file in altitude whereas the DSMC calculations do not include any diabatic 

849 heating. However, the results of Volkov et al. ( 2011a|[b ) also indicate that 

850 consistency with kinetic theory can be enforced approximately by applying 

851 the modified Jeans conditions to the hydrodynamic model at some altitude 

852 close to the exobase. This result is likely to be more general, and it applies 

853 to ionized atmospheres as long as ambipolar diffusion is taken into account 



854 (see Section 2.1.2). 

855 We compared the temperature and velocity profiles from the CI and C2 

856 models with results from similar models C5 and C6 that use the modified 

857 Jeans conditions. Note again that our version of the modified Jeans condi- 

858 tions includes the polarization electrostatic field that is required in strongly 

859 ionized media. Figure [8] shows the temperature and velocity profiles from the 

860 models. There is no difference between the C5 model and the CI model as 

861 long as the upper boundary of the C5 model is at a sufficiently high altitude. 

862 In this case we extended it to 36 Rp. When the upper boundary is placed at 

863 lower altitudes, the flow decelerates and the temperature increases near the 

864 upper boundary. A comparison between the C2 and C6 models provides an 
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865 example of the difference that can arise when the modified Jeans boundary 

866 conditions are used significantly below the exobase. A better agreement is 

867 achieved if the upper boundary of the C6 model is placed at a slightly higher 

868 altitude. In summary, we have shown that the CI and C2 models are both 

869 consistent with kinetic theory. 

870 We note that extending the models to 16 Rp or higher is not necessarily 

871 justified because it ignores the complications arising from the possible influ- 

872 ence of the stellar tide, the stellar wind, and interactions of the flow with the 

873 magnetic field of the planet or the star. We placed the upper boundary at 

874 a relatively high altitude to make sure that the boundary is well above the 

875 region of interest, but generally we do not consider our results to be accurate 

876 above 3-5 Rp. Instead, our results provide robust lower boundary conditions 

877 for multidimensional models of the escaping material outside the Roche lobe 

878 of the planet. Such models often cannot include detailed photochemical or 

879 thermal structure calculations. The results from the more complex models 

880 can then be used to constrain the upper boundary conditions in ID models. 

881 This type of an iteration is a complex undertaking, and it will be pursued in 

882 future work. 

883 3.2. Density profiles 

884 In this section we provide a qualitative understanding of the density pro- 
ses files and transition altitudes that affect the interpretation of the observations. 
886 Based on the gas giants in the solar system it might be expected that heavy 
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Figure 8: Temperature (a) and velocity (b) profiles in the upper 
atmosphere of HD209458b based on models with extrapolated 
and modified Jeans upper boundary conditions (see Table [2] for 
the input parameters). 
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species undergo diffusive separation in tlie tliermospfiere. If tliis were tlie 



888 case, tlie transit deptlis in tlie O I, C II, and Si III lines ( Vidal-Madjar et al. 



2004 Linsky et al. , 2010) should not be significantly higher than the transit 



890 depth at visible wavelengths. It is therefore important to explain why diffu- 

891 sive separation does not take place in the thermosphere of HD209458b, and 

892 to clarify why H and O remain mostly neutral while C and Si are mostly ion- 

893 ized. Also, doubly ionized species such as Si^"*" are not common in planetary 

894 ionospheres, and their presence needs to be explained. In order to do this, 

895 we modeled the ionization and photochemistry of the relevant species, and 

896 prove that diffusive separation does not take place. 

897 In order to illustrate the results. Figure |9] shows the density profiles of 

898 H, H+, He, He+, O, 0+, C, C+, Si, Si+, and Si2+ from the C2 model. The 

899 location of the H/H"*" transition obviously depends on photochemistry, but it 

900 also depends on the dynamics of escape. With a fixed pressure at the lower 

901 boundary, a faster velocity leads to a transition at a higher altitude. Thus 

902 the transition occurs near 3.1 Rp in the C2 model whereas in the CI and C3 

903 models it occurs at 3.8 Rp and 5 Rp, respectively. These results disagree with 



Yelle (2004) and Murray-Clay et al. (2009) who predicted a lower transition 



905 altitude, but they agree qualitatively with the solar composition model of 



Garcia Munoz (2007). They also agree with the empirical constraints derived 



907 by Koskinen et al. (2010a) from the observations. 



Once again, the differences between the earlier models and our work arise 
from different boundary conditions, and assumptions regarding heating rates 
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Figure 9: Density profiles in the upper atmosphere of 
HD209458b based on the C2 model (see Table g for the input 
parameters). 
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910 and photochemistry. We demonstrate this by reproducing the results of 



Murray-Clay et al. (2009) with our model. In order to do so, we set the 



lower boundary to 30 nbar with a temperature of 1,000 K, and included 
the substellar tide in the equations of motion. We only included H, H+, 
and electrons in the model, and used the recombination rate coefficient and 



Lyman a cooling rate from Murray-Clay et al. (2009). We also calculated 



the heating and ionization rates with the gray approximation by assuming a 
single photon energy of 20 eV for the stellar flux of 0.45 W m~^ at the orbital 



position of HD209458b. Figure 10 shows the density profiles of H and H"*" 



based on this model (hereafter, the MC09 model). 

The H/H+ transition in the MC09 model occurs near 1.4 Rp. If we re- 
place the gray approximation with the full solar spectrum in this model, the 
II/II+ transition moves higher to 2-3 Rp. This is because photons with dif- 
ferent energies penetrate to different depths in the atmosphere, extending 
the heating profile in altitude around the heating peak. This is why the tem- 
perature at the 30 nbar level in the C2 model is 3,800 K and not 1,000 K. 
In order to test the effect of higher temperatures in the lower thermosphere, 
we extended the MC09 model to = 1 yubar (with Tq = 1,300 K) and again 
used the full solar spectrum for heating and ionization. With these condi- 
tions, the H/II+ transition moves up to 3.4 Rp, in agreement with the C2 
model. We conclude that the unrealistic boundary conditions and the gray 



931 approximation adopted by Murray-Clay et al. (2009) and Guo (2011) lead 



932 to an underestimated overall density of H and an overestimated ion fraction. 
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933 Thus their density profiles yield a H Lyman a transit depth of the order of 

934 2-3 % i.e., not significantly higher than the visible transit depth. 



935 We note that Yelle (2004) also predicted a relatively low altitude of 1.7 

936 R.p for the H/H"*" transition - despite the fact that his model does not rely on 

937 the gray approximation and the lower boundary is in the deep atmosphere. 

938 The reason for the low altitude of the H/H^ transition in this case is the 

939 neglect of heavy elements. In the absence of heavy elements, H;^ forms near 

940 the base of the model and subsequent infrared cooling balances the EUV 

941 heating rates. This prevents the dissociation of H2 below the 10 nbar level. 

942 In reality, reactions with OH and thermal decomposition dissociate H2 near 



943 the 1 yubar level (see Section 2.1.1) and cooling by Kg is negligible at all 

944 altitudes. It should be noted that even if II2 does not initially dissociate, Hj^ 

945 can be removed from the lower thermosphere in reactions with carbon and 



946 oxygen species (e.g., Garcia Munoz, 2007) unless these species undergo diffu- 

947 sive separation. The subsequent lack of radiative cooling will then dissociate 

948 II2 again near the 1 /ibar level. 

949 In our models, charge exchange with oxygen (reactions R14 and R15 in 

950 Table [1]) dominates the photochemistry of H below 3 Rp and charge exchange 

951 with silicon (R25, R26) is also important below 1.4 Rp. These reactions are 

952 secondary in a sense that they require the ions to be produced by some other 

953 mechanism. In fact, 11+ is mainly produced by photoionization (PI), al- 

954 though thermal ionization (R3) is also important near the temperature peak. 

955 The production rates are mainly balanced by loss to radiative recombination 
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Figure 10: Density profiles of H and based on tlie MC09 



model that is similar to that of Murray-Clay et al. (2009) (see 
text). Compared with our models (see Figure 9), the H/H+ 
transition occurs at a significantly lower altitude. The difference 
arises from the lower boundary conditions and gray approxima- 
tion to heating and ionization in the MC09 model. 
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956 (Rl). The net chemical loss timescale for H is longer than the timescale for 

957 advection above 1.7 Rp. This allows advection from below to replenish H at 

958 higher altitudes. 

959 The density profiles of O and 0+ are strongly coupled to H and H"*" 



960 by charge exchange (see also Garcia Munoz, 2007). As a result, the 0/0"*" 

961 transition occurs generally near the H/H^ transition. For instance, in the 

962 C2 model it is located near 3.4 Rp. The same is not true of carbon. The 

963 C/C^ transition occurs at a much lower altitude than the H/H"*" and 0/0^ 

964 transitions. For instance, in the C2 model it is located near 1.2 Rp. is 

965 mainly produced by photoionization (P4), although thermal ionization (R8) 

966 and charge exchange with He+ (R13) are also important near the temperature 

967 peak. The production is balanced by loss to radiative recombination (RIO). 

968 The chemical loss timescale for C is shorter than the timescale for advection 

969 below 1.8 Rp. Thus advection is unable to move the C/C+ transition to 

970 altitudes higher than 1.2 Rp. 

971 Silicon is almost fully ionized near the lower boundary of the model. Much 

972 of the Si"*" below 4 Rp is produced by charge exchange of Si with H"*", He"*", and 

973 (R22, R23, R24). The low ionization potential of Si (8.2 eV) means that 

974 Si^ can also be produced by thermal ionization (R18), and photoionization 

975 (P6) by stellar FUV radiation and X-rays that propagate past the EUV 



976 ijieating peak. Above 4 Rp, Si"^ is mostly produced by photoionization. Linsky 



et al. (2010) suggested that the balance of Si"*" and Si^"*" depends on charge 



978 exchange with H"*" and H, respectively, and our results confirm this. However, 
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979 the location of the Si"'"/Si^"'" transition also depends on the dynamics. For 

980 instance, in the C2 model it occurs near 5.8 Rp while in the CI model it 

981 occurs near 8.5 Rp. Thus slow outflow and high temperatures favor Si^+ as 

982 the dominant silicon species as long as the flux constant is high enough to 

983 overcome diffusive separation (Paper II). 

984 We have now explained the presence of the atoms and ions that have been 

985 detected in the thermosphere of HD209458b. Due to advection and charge 

986 exchange, H and O are predominantly neutral up to about 3 -Rp and give rise 

987 to the observed transit depths in the H Lyman a and O I lines. Carbon, on 

988 the other hand, is ionized at a low altitude and thus C"*" is also detectable in 

989 the upper atmosphere. Si"^ is the dominant sihcon species below 5 Rp, but 

990 charge exchange with H ensures that there is also a significant abundance of 

991 Si^^ in the atmosphere. We note that these conclusions are only valid if the 

992 heavier species are carried along to high altitudes by the escaping hydrogen. 

993 We show that this is the case below in Section 13.2.21 

994 3.2.1. The EUV ionization peak (EIP) layer 



995 Koskinen et al. (2010b) explored the properties of the ionospheres of 

996 EGPs at different orbital distances from a Sun-like host star by using a hy- 

997 drostatic general circulation model (GCM) that also includes realistic heat- 

998 ing rates, photochemistry, and transport of constituents. They predicted 

999 that the EIP layer on HD209458b is centered at 1.35 Rp where the electron 

1000 density is rif. = 6.4 x 10^'^ m~^ and ~ 3 x 10^^. In the C2 model, the 
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EIP layer is centered at 1.3 Kp {p = 2 nbar) where = 4.4 x 10^^ m ^ and 



1002 Xp 



3.7 X 10 . The vertical outflow velocity at 1.3 R„ is 90 m s . Thus 



1003 the results of Koskinen et al. (2010b) were not significantly affected by the 

1004 simplifying assumptions of the GCM. This means that hydrostatic GCMs 

1005 can be extended to relatively low pressures as long as the escape rates are 

1006 incorporated as boundary conditions. 

1007 We also calculated the plasma frequency based on the electron densities 

1008 in the C2 model. This constrains the propagation of possible radio emissions 

1009 from the ionosphere. The ordinary plasma frequency Up/2Tr exceeds 12 MHz 

1010 at all altitudes below 5 Rp and reaches a maximum of almost 64 MHz in 

1011 the EIP layer. This presents a limitation on the detection of radio emissions 

1012 from the ionospheres of close-in EGPs. Any emissions that originate from 

1013 the ionosphere at 1-5 Rp and have frequencies lower than 10-70 MHz can 

1014 be screened out by the ionosphere itself. We note that a detection of radio 



1015 emissions from an EGP has not yet been achieved (e.g., Bastian et al. , 2000 



Lazio et al. 2007; Lecavelier des Etangs et al. 2011 ; Griefimeier et al. 2011 ). 



Such a detection would be an important constraint on the magnetic field 



1018 strength and the ionization state of the source region (e.g., Griefimeier et al. 



2007). Models of the ionosphere are required to predict radio emissions from 



the possible targets. 
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1021 3.2.2. The escape of heavy atoms and ions 

1022 In this section we verify a posteriori that the velocity and temperature 

1023 differences between different species in the thermosphere of HD209458b are 

1024 smalL This demonstrates that the single fluid approximation of the mo- 

1025 mentum and energy equations is valid, and that diffusive separation of the 

1026 heavy species does not take place. Our model includes velocity differences 

1027 between different species by including all of the relevant collisions between 

1028 them through the inclusion of diffusive fluxes in the continuity equations. 

1029 However, we have explicitly assumed that Tn — Ti — Tg, and this assumption 

1030 in particular needs to be verified. Also, the diffusion approach to the con- 

1031 tinuity equation is only valid if the velocity differences between the species 

1032 are reasonably small. 

1033 We calculated the collision frequencies based on the C2 model, and found 

1034 that collisions with neutral H dominate the transport of heavy neutral atoms 

1035 such as O below 3.5 Rp. At altitudes higher than this, collisions with H+ are 

1036 more frequent. In Paper II we demonstrate that a mass loss rate of 6 x 10^ 

1037 kg s~^ is required to prevent diffusive separation of O (the heaviest neutral 

1038 species detected so far) in the thermosphere. The mass loss rate in our mod- 

1039 els is M > 10^ kg s~^ and thus O is dragged along to high altitudes by H. On 

1040 the other hand, collisions with H"*" dominate the transport of heavy ions such 

1041 as Si"*" as long as the ratio [H+]/[H] > 10~^ (Paper II). This explains why 

1042 Coulomb collisions in our models are more frequent than heavy ion-H coUi- 

1043 sions at almost all altitudes apart from the immediate vicinity of the lower 
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1044 boundary. These collisions are much more efficient in preventing diffusive 

1045 separation than collisions with neutral H. 



1046 Figure 11 compares the timescale for diffusion td for O and Si+ with the 

1047 timescale for advection based on the C2 model. In both cases, td » 

1048 and thus diffusion is not significant. This implies that there are no signifi- 

1049 cant velocity differences between heavy atoms and hydrogen. We note that 

1050 Coulomb collisions of doubly ionized species with are more frequent than 

1051 collisions between a singly ionized species and H+. Thus diffusion is even 

1052 less significant for a species like Si^"^. We verified these results from our sim- 

1053 ulations by switching diffusion off in the model and rerunning the C2 model. 

1054 As a result the density of the heavy atoms and ions increased slightly at high 

1055 altitudes, but the differences are not significant - the results were nearly 

1056 identical to the density profiles of the original simulation. 

1057 We note that the atmosphere can also be mixed by vertical motions asso- 

1058 ciated with circulation that are sometimes parameterized in one-dimensional 



1059 models by the eddy diffusion coefficient Kzz (e.g., Moses et al. , 2011). This 

1060 mechanism is efficient in bringing the heavy elements to the lower thermo- 

1061 sphere but it is unlikely to mix the atmosphere up to 3 Rp and beyond. Also, 

1062 there is no generally accepted method of estimating the degree of global mix- 

1063 ing based on circulation models, and most circulation models for EGPs do 

1064 not adequately treat the relevant energy deposition and forcing mechanisms 

1065 in the upper atmosphere. Thus there is considerable uncertainty over the val- 

1066 ues of Kzz and rapid escape is a much more likely explanation for the lack of 
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diffusion 
Si* diffusion 
Advection 



1 10 

Figure 11: Timescales for diffusion td = H^/Dg of O and 
Si"^, and for advection = H/v based on the C2 model. 
We calculated the diffusion coefficients in a mixture of H and 
H+. The large scale height of the atmosphere and relatively 
high collision frequencies mean that diffusion is not significant 
{td/tv = vH/Ds >> 1) in the thermosphere of HD209458b. 
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1067 ijiiffusive separation on HD209458b. In fact, the calculations of Koskinen et 



1068 |aL] ( |2007b ) show that the temperature in the thermosphere of planets such 



1069 as HD209458b is high enough to practically guarantee an effective escape 

1070 rate. The only way to prevent this is to provide enough radiative cooling to 

1071 offset most of the XUV flux, but there are no radiative cooling mechanisms 

1072 efficient enough to achieve this in a thermosphere composed of atoms and 

1073 ions. 

1074 As we noted above, the temperatures of the electrons, ions, and neutrals 

1075 are roughly equal in the thermosphere of HD209458b. In order to show this, 

1076 we assumed that photoelectrons share their energy with thermal electrons, 

1077 which then share this energy further with ions and neutrals. We also as- 

1078 sumed that the collisions frequencies between the species are higher than 

1079 the timescale for advection. If the velocity differences between the species 

1080 are negligible, the steady state 5-moment energy equations for thermal elec- 



1081 trons and ions (Schunk and Nagy, 2000) can be used to obtain the following 
approximations^ 



Te-Ti ^ l^i^ (15) 

3 rrie krieVei 

Ti-Tn ^ (16) 

3 Uli KTliUin 



where qr is the volume heating rate, and z/gj and are the electron-ion and 



''Note that conduction and viscosity are not important in the thermosphere of 
HD209458b. 
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ion-neutral momentum transfer collision frequencies, respectively. 

We calculated the temperature differences for H, H+, and electrons based 
on the C2 model. The difference between the electron and ion temperatures 
decreases with altitude and is mostly less than 2 K. The difference between 
the ion and neutral temperatures, on the other hand, increases with altitude. 
The ion temperature is approximately 10 K higher than the neutral temper- 
ature near 5 Rp and the difference reaches 150 K at 16 Rp. In both cases, the 
temperature differences are negligible compared to the temperature of the 
thermosphere. Further, the timescale for advection in the C2 model is al- 
ways significantly longer than the relevant collision timescales. Thus we have 



shown that Tg ~ T„ and that equations (15) and (16) are approximately 
valid. 



1096 4. Discussion and Conclusions 



1097 We have constructed a new model for the upper atmosphere of HD209458b 

1098 in order to explain the detections of H, O, C"*", and Si^"^ at high altitudes 



1099 around the planet ( Vidal-Madjar et al., 2003, 2004 Linsky et al., 2010) 



1100 There are many different interpretations of the observed transits in the H 



2003 


Ben- Jaffel 


2007, 


2008 


Holstrom 



et al. , 2008; Koskinen et al. , 2010a), and these interpretations rely on results 



1103 from models of the upper atmosphere that are based on many uncertain as- 



1104 sumptions (see Section 3.1.1 and Koskinen et al. , 2010a, for a review). Also 



1105 the detection of heavy atoms and ions in the thermosphere is not without 
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controversy, and the detection of Si^^ is particularly intriguing. Thus these 
observations present several interesting challenges to modelers. 

The observed transit depths are large, and substantial abundances of the 
relevant species are required to explain the observations. However, on every 
planet in the solar system heavier species are removed from the thermo- 
sphere by molecular diffusion and doubly ionized species are not commonly 
observed. Also, the observations imply that H and O remain mostly neutral 
in the thermosphere while C and Si are mostly ionized at a relatively low 
altitude. Hydrodynamic models coupled with chemistry and thermal struc- 
ture calculations are required to explain the detection of these species in the 
upper atmosphere and the differences between their density profiles. Ours is 
the first such model that benefits from repeated detections of both neutral 
atoms and ions to constrain the composition and temperature. 



Koskinen et al. (2010a) demonstrated that the H Lyman a transit obser- 



1120 vations (Ben-Jaffel, 2007, 2008) can be explained with absorption by H in 



the thermosphere if the base of the hot layer of H is near 1 /xbar, the mean 
temperature within the layer is about 8,250 K, and the atmosphere is mostly 
ionized above 3 R^. These parameters are based on fitting the data with 
a simple empirical model of the upper atmosphere. The density and tem- 
perature profiles from our new hydrodynamic model agree qualitatively with 



ifhese constraints, demonstrating that the basic assumptions of [Koskinen et 



1127 al. (2010a) are reasonable. This confirms once again that a comet-like tail 



1128 (Vidal-Madjar et al. 2003) or energetic neutral atoms (Holstrom et al. 2008) 
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1129 are not necessarily required to explain the H Lyman a observations. 



1130 In line with recent results by Moses et al. (2011) and the empirical con- 
mi straints mentioned above, we used a photochemical model of the lower atmo- 

1132 sphere to show that H2 dissociates near the 1 /ibar level. Above this level, 

1133 the lack of efficient radiative cooling and strong stellar EUV heating lead 

1134 to high temperatures. We constrained the range of possible mean (pressure 

1135 averaged) temperatures based on the average solar flux by using the hydro- 

1136 dynamic model to calculate temperatures with different heating efficiencies. 

1137 For net heating efficiencies between 0.1 and 1, the mean temperature below 

1138 3 -Rp varies from 6,000 K to 8,000 K. This means that 8,000 K is a relatively 

1139 strict upper limit on the mean temperature if the XUV flux of HD209458 is 

1140 similar to the corresponding flux of the sun. 

1141 A mean temperature of 8,250 K estimated from the observations implies 

1142 the presence of an additional non-radiative heat source, or that the XUV flux 

1143 from HD209458 is higher than the average solar flux. Given that our best 



1144 estimate of the net heating efficiency is 0.44 (see Section 3.1.2), the XUV flux 

1145 of H209458 would have to be 5-10 times higher than the average solar flux to 



1146 cause a mean temperature of 8,250 K (see Section 3.1.1). If the mean XUV 

1147 flux of HD209458 is generally higher than the solar flux and the observations 

1148 took place during stellar maximum, such an enhancement is not impossible. 

1149 This would also lead to higher outflow velocity and mass loss rate. However, 



1150 the uncertainty in the observed transit depths is also large (Ben-Jaffel, 2008 



Ben-Jaffel and Hosseini 2010), and it can accommodate a range of temper- 
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1152 atures. Therefore our reference model C2 with a mean temperature of 7,200 

1153 K also agrees qualitatively with the empirical constraints. In this respect, it 

1154 is interesting to note that with lOOx solar flux, the mean temperature is still 



1155 ipnly 9,800 K. Temperatures significantly higher than 8,000 K (e.g., Ben-Jaffel 



and Hosseini, 2010) therefore imply a strong non-radiative heat source. 



1157 In contrast to the mean temperature, the velocity and details of the tem- 

1158 perature profile depend strongly on the heating efficiency and stellar flux 



1159 (see Section 3.1.1). They are also sensitive to the upper and lower boundary 

1160 conditions. This explains the large range of temperature and density profiles 

1161 predicted by earlier models that arise from different boundary conditions and 

1162 assumptions about the stellar flux, radiative transfer, and heating efficiencies. 

1163 The differences highlight the need for accurate thermal structure calculations 

1164 that are constrained by the available observations. These calculations are im- 

1165 portant because the density profiles of the detected species depend on the 

1166 temperature and velocity profiles, and inappropriate assumptions made by 

1167 the models can bias the interpretation of the observations. 

1168 In the absence of stellar gravity, the location of the sonic point and the 

1169 outflow speed also depend on the heating efficiency. As the heating efficiency 

1170 increases from 0.1 (in models with the average solar flux), the high altitude 

1171 temperature increases and the sonic point moves to lower altitudes, reaching 

1172 down to 4 Rp with a net heating efficiency of 1. We found that supersonic 

1173 solutions are possible as long as there is significant heating over a large alti- 

1174 tude range above the temperature peak. This conclusion is supported both 
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1175 |j3y the hydrodynamic model and new constraints from kinetic theory ( Volkov 



et aL, 2011a|[b ). However, the isentropic sonic point of the C2 model is lo- 



1177 cated above the model domain. In principle, this is an interesting result but 

1178 it should be treated with caution. We used parameterized heating efficiencies 

1179 for low energy photons, and the location of the sonic point is very sensitive to 

1180 the temperature profile. Also, the stellar tide can enhance the escape rates at 

1181 the substellar and antistellar points. We did not include this effect because 

1182 it may produce horizontal flows that cannot be modeled in ID. 

1183 As long as the upper boundary is at a sufficiently high altitude, we found 

1184 that the results based on the outflow boundary conditions and modifled Jeans 



1185 conditions are identical (see Section 3.1.5). This shows that our simulations 

1186 are roughly consistent with kinetic theory. An agreement between these two 

1187 types of boundary conditions on HD209458b is an interesting theoretical 

1188 result. It shows that the boundary conditions for hydrodynamic escape are 

1189 appropriate in this case. However, an upper boundary at 16 Rp or higher 

1190 is not necessarily justifled for other reasons because we did not consider the 

1191 effect of the possible planetary magnetic fleld, interaction of the atmosphere 



1192 with the stellar wind, or horizontal transport (e.g.. Stone and Proga, 2009 



Trammell et al. , 2011). 



1194 We chose an upper boundary at a high altitude in order to preserve the 

1195 integrity of the solution in our region of interest below 5 Rp. The purpose of 

1196 this work is to model energy deposition and photochemistry in this region. 

1197 These aspects are often simplifled in more complex models to a degree that 
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1198 it may be difficult to separate the effect of multiple dimensions and other 

1199 complications from differences arising simply because of different assump- 

1200 tions about heating efficiencies and chemistry. Also, the uncertainty in the 

1201 observations does not necessarily justify the introduction of more free pa- 

1202 rameters to the problem until the basic properties of the thermosphere are 

1203 better understood. However, technically we do not consider our solutions 

1204 to be accurate far above 3-5 Rp. Instead, our results provide robust lower 

1205 boundary conditions for more complex multidimensional models that char- 

1206 acterize the atmosphere outside the Roche lobe of the planet. Results from 

1207 such models can then be used to constrain the upper boundary conditions of 

1208 the ID models further. 

1209 In order to model the density proffies of the detected species in the iono- 



1210 sphere, we assumed solar abundances of the heavy elements (Lodders, 2003), 

1211 although this assumption can be adjusted as required to explain the obser- 

1212 vations (Paper II). As we already stated we found that H2, H2O, and CO 

1213 dissociate above the 1 /xbar level, releasing H, O, and C to the thermosphere 



1214 (see also Moses et al. , 2011 ). We note that the detection of Si^+ in the upper 

1215 atmosphere implies that silicon does not condense into clouds of forsterite 



1216 and enstatite in the lower atmosphere as argued by e.g., Visscher et al. (2010 ). 

1217 The dominant Si species is then SiO, which dissociates at a similar pressure 

1218 level as the other molecules. In fact, practically all molecules dissociate below 

1219 0.1 /ibar. This leads to an important simplification in hydrodynamic models 

1220 of the thermosphere. The complex chemistry of molecular ions does not need 
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1221 to be included as long as the lower boundary is above the dissociation level. 

1222 We found that the H/H+ transition occurs near 3 Rp or, depending on 

1223 the velocity profile, at even higher altitudes. The 0/0"*" transition is coupled 

1224 to the H/H+ transition through charge exchange reactions. Thus both H and 

1225 O are mostly neutral up to the boundary of the Roche lobe at 3-5 Rp. In 

1226 contrast, C is ionized near the 1 /xbar level and C"*" is the dominant carbon 

1227 species in the thermosphere. Si is also ionized near the 1 //bar level, and the 

1228 balance between Si"*" and Si^"*" is determined by charge exchange with H"*" and 

1229 H, respectively. Si"*" is the dominant silicon ion below 5 Rp but the abundance 

1230 of Si^+ is also significant. We found that neutral heavy atoms are dragged 

1231 to the thermosphere by the escaping H, while heavy ions are transported 

1232 efficiently by the escaping H+. Thus the advection timescale is much shorter 

1233 than the diffusion timescale of the detected species, and diffusive separation 

1234 does not take place in the thermosphere. We also verified that the neutral, 

1235 ion, and electron temperatures are roughly equal. 

1236 Taken together, these results imply that the thermospheres of close-in 

1237 EGPs can differ fundamentally from the gas giant planets in the solar sys- 

1238 tem. For instance, the thermosphere of HD209458b is composed mainly of 

1239 atoms and atomic ions, and diffusive separation of the common heavy species 

1240 is prevented by the escape of H and H+. It is important to note, however, 

1241 that results such as these cannot be freely generalized to other extrasolar 

1242 planets. As in the solar system, each planet should be studied separately. 

1243 For instance, the dissociation of molecules depends on the temperature pro- 
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1244 file that is shaped by the composition through radiative coohng and stellar 

1245 heating. The mass loss rate and escape velocity, that determine whether dif- 

1246 fusive separation takes place or not, depends on the escape mechanism that 

1247 again depends on the temperature and composition of the upper atmosphere. 

1248 The results from different models can only be verified by observations that 

1249 are therefore required for multiple planets if we are to characterize escape in 

1250 different systems and under different conditions. 
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